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ABSTRACT 

Origins of spectral methods, especially their relation to the Method of 
Weighted Residuals, are surveyed. Basic Fourier, Chebyshev, and Legendre 
spectral concepts are reviewed, and demonstrated through application to simple 
model problems. Both collocation and tau methods are considered. These 
techniques are then applied to a number of difficult, nonlinear problems of 
hyperbolic, parabolic, elliptic, and mixed type. Fluid-dynamical applications 
are empha sized. 
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INTRODUCTION 


Spectral methods may be viewed as an extreme development of the class of 
discretization schemes known by the generic name of the method of weighted 
residuals (MWR) [1]. The key elements of the MWR are the trial functions 
(also called the expansion or approximating functions) and the test functions 
(also known as weight functions)* The trial functions are used as the basis 
functions for a truncated series expansion of the solution, which, when 
substituted into the differential equation, produces the residual. The test 
functions are used to enforce the minimization of the residual. 

The choice of trial functions is what distinguishes the spectral methods 
from the finite element and finite difference methods. The trial functions 
for spectral methods are infinitely differentiable global functions. 
(Typically they are tensor products of the eigenfunctions of singular Sturm- 
Liouville problems.) In the case of finite element methods, the domain is 

divided into small elements, and a trial function is specified in each 
element. The trial functions are thus local in character, and well-suited for 
handling complex geometries. The finite difference trial functions are 
likewise local. 

The choice of test function distinguishes between the Galerkin, 
collocation, and tau approaches. In the Galerkin approach, the test functions 
are the same as the trial functions, whereas in the collocation approach the 
test functions are translated Dirac delta functions. In other words, the 

Galerkin approach is equivalent to a least squares approximation, whereas the 
collocation approach requires the differential equation to be satisfied 
exactly at the collocation points. Spectral tau methods are close to Galerkin 
methods but they differ in the treatment of boundary conditions. 
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The collocation approach is the simplest of the MWR, and appears to have 
been first used by Slater [2] in his study of electronic energy bands in 
metals. A few years later, Barta [3] applied this method to the problem of 
the torsion of a square prism. Frazer, et al. [4] developed it as a general 
method for solving ordinary differential equations. They used a variety of 
trial functions and an arbitrary distribution of collocation points. The work 
of Lanczos [5] established for the first time that a proper choice of trial 
functions and distribution of collocation points is crucial to the accuracy of 
the solution. Perhaps he should be credited with laying down the foundation 
of the orthogonal collocation method. This method was revived by Clenshaw [6] , 
Clenshaw and Norton [7], and Wright [8]. These studies involved application 
of Chebyshev polynomial expansions to initial value problems. Villadsen and 
Stewart [9] developed this method for boundary value problems. 

The earliest investigations of the spectral collocation method to partial 
differential equations were those of Kreiss and Oliger [10] (who called it the 
Fourier method) and Orszag [11] (who termed it pseudospectral) . This approach 
is especially attractive because of the ease with which it can be applied to 
variable coefficient and even nonlinear problems. The essential details will 
be furnished below. 

The Galerkin approach is perhaps the most esthetically pleasing of the 
MWR since the trial functions arid the test functions are the same. Indeed, 
the first serious application of spectral methods to FDE s — — that of 
Silberman [12] for meteorological modelling — used the Galerkin approach. 
However, spectral Galerkin methods only became practical for high resolution 
calculations of nonlinear problems after Orszag [13] and Eliasen, et al. [14] 
developed a transform method for evaluating convolution sums arising from 
quadratic nonlinearities. Even in this case spectral collocation methods 
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retain a factor of 2 in speed* For more complicated nonlinear terms high 
resolution spectral Galerkin methods are still impractical* 

The tau approach is the most difficult to rationalize within the context 
of the MWR. Lanczos [5] developed the spectral tau method as a modification 
of the Galerkin method for problems with non-periodic boundary conditions. 
Although it too, is difficult to apply to nonlinear problems, it has proven 
quite useful for constant coefficient problems or subproblems, e.g*, for semi- 
implicit time-stepping algorithms* 

The following discussion of spectral methods for FDE's will be organized 
around the three basic types of systems — hyperbolic, parabolic, and elliptic 
— with an additional section for a difficult, nonlinear problem of mixed 
type* Simple, one-dimensional, linear examples will be provided to illustrate 
the basic principles and details of the algorithms; two-dimensional, nonlinear 
examples drawn from fluid dynamical applications will also be furnished to 
demonstrate the power of the method. The focus will be on collocation 
methods, although some discussion of tau methods is provided. 


II. HYPERBOLIC EQUATIONS 

Linear hyperbolic equations are perhaps the simplest setting for 
describing spectral collocation methods* Both Fourier and Chebyshev schemes 
have found wide application. This section will first present the fundamentals 
of both approaches and then illustrate them on a nonlinear fluid dynamics 
problem involving shock waves. 
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Basic Fourier Spectral Concepts . 

The potential accuracy of spectral methods derives from their use of 
suitable high-order interpolation formulae for approximating derivatives. An 
elementary example is provided by the model problem 

u t + u x = 0, (1) 

with periodic boundary conditions on [0 ,2 tt ] and the initial condition 

u(x,0) = sin(ir cos x). (2) 

The exact solution 

u(x,t) = sin[TT cos(x-t)] ( 3 ) 

has the Fourier expansion 

u(x,t) = l u k (t) e ikX , (M 

k=-°° 

where the Fourier coefficients 

\(t) = sin(^) J k (7T) e" ikt (5) 

and J k (t) is the Bessel function of order k. The asymptotic properties of 
the Bessel functions imply that 

k p u k (t) - 0 


as k + 00 


( 6 ) 
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for all positive integers p» As a result, the truncated Fourier series 


u N (x,t) 


N/2-1 

l 

k- -N/2+1 


u k ( t) e 


ikx 


(7) 


converges faster than any finite power of 1/N. This property is often 
referred to as exponential convergence. A straightforward integration— by- 
parts argument [15] may be used to show that it applies to any periodic and 
infinitely differentiable solution. 

The standard collocation points are 


x . 
J 



j=0, l,***, N-l. 


( 8 ) 


Let uj denote the approximation to u(xj), where the time dependence has 
been suppressed. Then the spatial discretization of Eq. (1) is 




(9) 


where the right— hand— side is determined as follows. First, compute the 
discrete Fourier coefficients 


N- 

I u 

J-o 


-ikx . 
e J 


Then the interpolating function 




N 

2 


1 . 


( 10 ) 


u(x) 


N/2-1 


l u,. e 


ikx 


k= -N/2 


( 11 ) 


can be differentiated analytically to obtain 
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111 
9x j 


N/2-1 

I 

k» -N/2+1 


ik °k 


ikx . 
e J 


( 12 ) 


(The term involving k m — N/2 makes a purely imaginary contribution to the 
sum and hence has been dropped.) Note that each derivative approximation uses 
all available information about the function values. The sums in Eqs. (10) 
and (12) can be obtained in 0(N in N) operations via the Fast Fourier 
Transform (FFT). 

An illustration of the superior accuracy available from the spectral 
method for this problem is provided in Table I. Shown there are the maximum 
errors at t = 1 for the truncated series and for the spectral collocation 
method as well as for second-order and fourth-order finite difference 
methods. The time discretization was the classical fourth-order Runge-Kutta 
method. In all cases the time-step was chosen so small that the temporal 
discretization error was negligible. Because the solution is infinitely 
smooth, the convergence of the spectral method on this problem is more rapid 
than any finite power of 1/N. (The error for the N = 64 spectral result is 


Table I. Maximum Error for a 1-D Periodic Problem 


N 


2nd-0rder 4th-0rder 

Truncated Fourier Finite Finite 

Series Spectral Difference Difference 


8 

9.87 

(-2) 

16 

2.55 

(-4) 

32 

1.05 

(-11) 

64 

6.22 

(-13) 

128 




1.62 (-1) 1.11 

4.97 (-4) 6.13 

1.03 (-11) 1.99 

9.55 (-12) 5.42 

1.37 


(0) 

9.62 

(-1) 

(-1) 

2.36 

(-D 

(-1) 

2.67 

(-2) 

(-2) 

1.85 

(-3) 

(-2) 

1.18 

(-4) 
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so small that it is swamped by the round-off error of these single precision 
CDC Cyber 175 calculations.) In most practical applications the benefit of 
the spectral method is not the extraordinary accuracy available for large N 
but rather the small size of N necessary for a moderately accurate solution. 


Basic Chebyshev Spectral Concepts 

Spectral methods for non-periodic problems can also exhibit exponential 
convergence. A simple example is again provided by Eq. (1) but now on the 
interval [-1,1] with initial condition u(x,0) and boundary condition 

Since this is not a periodic problem, a spectral method based upon 
Fourier series in x would exhibit extremely slow convergence. However, 
rapid convergence as well as efficient algorithms can be attained for spectral 
methods based upon Chebyshev polynomials. These are defined on [-1,1] by 

T n (x) = cos (n cos *x). ( 13 ) 


The function 


u (x , t) = sin onr(x-t) 


(14) 


is one solution to Eq. (1). It has the Chebyshev expansion 


u(x,t) = l u (t) t (x), 
n=0 n n 


(15) 


where 

u n (t) sin (• s |- art) J n (ar) 

n 


with 


(16) 
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n * 0 


(17) 


c =* 
n 


i 2 

i. 


n > 1 


The truncated series 


u N (x,t) 


Jo “" (t) T " (x) 


(18) 


converges at an exponential rate. Note that this result holds whether or 
not a is an integer. In contrast, the Fourier coefficients of u(x,t) are 

— _ 1 „ictnt sin ir(a+k) __ i -icmt sin ir(a-k) 

V t; 2? 6 S+k 27 6 ^k * (19) 


For non-integer a these decay extremely slowly. 
The change of variables 


x = cos 0, 


( 20 ) 


the definition 

v(9,t) = u(cos e,t), (21) 

and Eq. (13) reduce Eq. (15) to 

oo 

v(0,t) = l u (t) cos nO. (22) 

n=0 

Thus, the Chebyshev coefficients of u(x,t) are precisely the Fourier 
coefficients of v(0,t)- This new function is automatically periodic- If 
u ( x , t ) is infinitely differentiable (in x) , then v(0,t) will be infinitely 
differentiable (in 0). Hence, straightforward integration-by-parts arguments 
lead to the conclusion that the Chebyshev coefficients of an infinitely 
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differentiable function will decay exponentially fast. Note that this holds 
regardless of the boundary conditions. 

A Chebyshev spectral method makes use of the interpolating function 

N 

u(x) - I u n T n ( x ). (23) 

n=0 


The standard collocation points are 


x. *= cos 


II 

N 


j = 0, N. 


(24) 


Thus, 


u - I co» T 1 . 

n=0 


(25) 


where Uj is the approximation to u(x^). The inverse relation is 


N 


l 


Nc n jio ' S 


u. COS • n ,1T ' ^ 


N 


n = 0, 1, 


(26) 


where 


j = 0 or N 
1 < j < N-l 


(27) 


The analytic derivative of this function is 


3u 

3x" 



T n (x), 


(28) 


where 


10 


-(1) 

u 

N+l 

-U) 

u 

N 


_ -d) 

c n u n 


0 


0 


a(1) 

Un +2 + 2 ( n+1 ) u n +l * 


n = N— 1 ,N— 2 , * * * , 0. 


( 29 ) 


(See [15] for the derivation of this recursion relation.) The Qiebyshev 
spectral derivatives at the collocation points are 

I . I cos — . (30) 

3 * 'a n-0 N 

Special versions of the FFT may be used for evaluating the sums in Eqs. (26) 
and (30). The total cost for a Chebyshev spectral derivative is thus 
0(N Jin N). 

The time— stepping scheme for Eq. (1) must use the boundary conditions to 
update u N (at x = -1) and the approximate derivatives from Eq. (30) to 
update uj for j-0, l t ,,# ,N-l. Note that no special formula is required for 
the derivative at j =0 (or x = +1). 

Results pertaining to a = 2.5 at t = 1 for a truncated Chebyshev 
series, a Chebyshev spectral method, a Fourier spectral method, and a second- 
order finite difference method are given in Table II. For this non-periodic 
problem Fourier spectral methods are quite inappropriate, but the Chebyshev 
spectral method is far superior to the finite difference method. 

The Chebyshev collocation points are the extreme points of t n (x). Note 
that they are not evenly distributed in x, but rather are clustered near the 
endpoints. The smallest mesh size scales as 1/N^« While this distribution 
contributes to the quality of the Chebyshev approximation and permits the use 
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of the FFT in evaluating the series, it also places a severe time-step 
limitation on explicit methods for evolution equations. 


Table II. Maximum Error for a 1-D Dirichlet Problem 


N 

Truncated 

Series 

Chebyshev 

Spectral 

Fourier 

Spectral 

Finite 

Difference 

4 

1.24 

(0) 

1.49 

(0) 

1.85 (0) 

1.64 

(0) 

8 

1.25 

(-D 

6.92 

(-1) 

1.92 (0) 

1.73 

(0) 

16 

7.03 

(-6) 

1.50 

(-4) 

2.27 (0) 

1.23 

(0) 

32 

1.62 

(-13) 

3.45 

(-11) 

2.28 (0) 

3.34 

(-1) 

64 

1.79 

(-13) 

9.55 

(-11) 

2.27 (0) 

8.44 

(-2) 


Application to Two-dimensional, Supersonic Flow 

Spectral methods have recently been applied successfully to the nonlinear 
hyperbolic system of equations which describes a two-dimensional inviscid gas 

[16.17] . The most serious complication over the simple model problems 

discussed above occurs when shock waves are present. If the shock occurs in 
the interior of the domain, then the truncated series for the discontinuous 
flow variables converges very slowly. Elaborate filtering strategies appear 
necessary to extract useful information from a calculation of such a situation 

[17.18] . This difficulty disappears, however, when the shock occurs at the 
boundary of the domain, as in shock-fitting as opposed to shock-capturing 
calculations. 

A schematic of the type of spectral shock-fitted calculations described 
below is illustrated in Fig. 1. At time t = 0 an infinite, normal shock 
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at x ■* 0 separates a rapidly moving, uniform fluid on the left from the 
fluid on the right which is in a quiescent state except for some specified 
fluctuation. The initial conditions are chosen so that in the absence of any 
fluctuation the shock moves uniformly in the positive x-direction with a Mach 
number (relative to the fluid on the right) denoted by Mg. In the presence 
of fluctuations the shock front will develop ripples. The shape of the shock 
is described by the function x g (y,t). The numerical calculations are used to 
determine the state of the fluid in the region between the shock .front and 
some suitable left boundary x^(t) and also to determine the motion and shape 
of the shock front itself. 

'Figure 1 is taken from a shock/ turbulence calculation [19] in which the 
downstream fluctuation is a plane vorticity wave that is periodic in y with 
period y^. Because of the initial value nature of the calculation, the fluid 
motion behind the shock is not periodic in x, as Fig. 1 makes abundantly 
clear. The interesting physical domain is given by 

x L (t) < x < x g (y,t) 

0 < y < y £ (31) 

t > 0. 


The change of variables 
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x - x L (t) 

x ~ x g (y,t) - *^(o 

Y = y/y £ (32) 

T = t 

produces the computational domain 

0 < X < 1 

0 < Y < 1 (33) 

T > 0. 

The fluid motion is modeled by the two-dimensional Euler equations* In 
terms of the computational coordinates these are 

Q x + B Q x + c Q Y = °» (34) 

YX 0 

y 

0 0 

U 0 

0 u 





and 




The contravariant velocity components are given by 

U = X + uX + vX 
t x y 

and (37) 

V - Y + uY + vY . 
t x y 

A subscript denotes partial differentiation with respect to the indicated 
variable. P, a, and S are all normalized by reference conditions at 

downstream infinity; u and v are velocity components in the x and y 

directions, both scaled by the characteristic velocity defined by the square 
root of the pressure-density ratio at downstream infinity. A value y = 1.4 
has been used. 

Let n denote the time level and At the time increment. The time 
discretization of Eq. (34) is 

Q = [ 1 — AtL n ] Q n (38) 

Q n+1 - y [Q n + (1 " AtL)Q] , (39) 

where L denotes the spatial discretization of The solution 


Q has the Chebyshev - Fourier series expansion 
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Q(X,Y,T) 


M 

l 


P“0 


N/2-1 

l 

q=-N/2 



(T) T p (C)e 


2*iqY 

» 


(40) 


where £ = 2X-1. The derivatives Q x and Qy are approximated by 


Q x - 2 I T Q^WW' 1 ' 1 . 
X p=0 q =-N/2 pq P 


(41) 


0,-2. ! f' 1 Q<»’ I) (T,x (5 )e 2 ’^, 

Y p=0 q=-N/2 pq P 


(42) 


where Q is computed from Q_._ in a manner analogous to Eq. (29), and 

pq r pq 


(M> 

pq 



(43) 


As a general rule the correct numerical boundary conditions for a 
spectral method are the same as the correct analytical boundary conditions* 
The global nature of the approximation avoids the need for special 
differentiation formulae at boundaries* At the same time spectral methods are 
quite unforgiving of incorrect boundary conditions. The inherent dissipation 
of these methods is so low that boundary errors quickly contaminate the entire 
solution. In many fluid dynamical applications the computational region must 
be terminated at some finite, artificial boundary. The difficulty at 
"artificial" boundaries is that analytically correct, fully nonlinear boundary 
conditions for systems are seldom known. One example of a workable artificial 
boundary condition for. the Euler equations is given in Ref. [20]. 

The most critical part of the calculation is the treatment of the shock 
front. The shock-fitting approach used here is desirable because it avoids 
the severe post-shock oscillations that plague shock-capturing methods. The 
time derivative of the Rankine-Hugoniot relations provides an equation for the 
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shock acceleration. This equation is integrated to update the shock position 
(see [20] for details). This method is a generalization of the finite 
difference method developed by Pao and Salas [21] for their study of the 
shock/vortex interaction. 

The nonlinear interaction of plane waves with shocks was examined at 
length in [19]. The numerical method used there was similar to the one 
described above but employed second-order finite differences in place of the 
present Chebyshev-Fourier spectral discretization. Detailed comparisons were 
made in [19] with the predictions of linear theory [22]. The linear results 
turned out to be surprisingly robust, remaining valid at very low (but still 
supersonic) Mach numbers and at very high incident wave amplitudes. The only 
substantial disagreement occurred for incident waves whose wave fronts were 
nearly perpendicular to the shock front. Hiis type of shock— turbulence 
interaction is a useful test of the spectral technique because the method can 

be calibrated in the regions for which linear theory has been shown to be 
valid . 

The most reliable numerical results can be obtained for the acoustic 
responses to acoustic waves. Unlike the vorticity responses, these require no 
differentiation of the flow variables, thus eliminating one extra source of 
error. Moreover, the acoustic reponse stretches much further behind the shock 
than the vorticity response, thus providing greater statistical reliability. 
Vorticity response results are reported in [23] . The incident pressure wave 
is taken to be 

i(kj*2. - “^t) 

Pi " A i e (44) 

where J^i = (ki^jk^y)', “ ® M g ai k^ x + ai ki and A£ is the amplitude. In 
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terms of the incidence angle Xi * (^1 cos ®i»^i s ^ n ® i ) • The linearized 
transmitted acoustic wave can be expressed in the same manner with all 
subscripts changed from 1 to 2. The amplification coefficient for the 
transmitted acoustic wave is then the ratio 

A'/A' . (45) 


Figure 2 indicates the transmission coefficient extracted from the 
computation. At each fixed value of X we perform a Fourier analysis in Y 
of the pressure. The Fourier coefficient for q = 1 provides the amplitude 
A^. In order to reduce the transients that would accompany an abrupt start of 
the calculation at full wave amplitude, an extra factor of s(t) is inserted 
into Eq. (44), where 


s(t) 


3(t/t g ) 2 - 2(t/t g ) 3 
1 


0 < t < t 

s 


t > t 

s 


(46) 


The start-up time t g is some multiple (typically V2 ) t ^ !e time takes 
the shock to encounter one full wavelength (in the x-direction) of the 
incident wave. The ratio A^/AJ is plotted in Fig. 2 as a function of the 
mean value of the physical coordinate x corresponding to X. The start-up 
time for this Mach 3 case is t g = 0.56. The average of the x-dependent 
responses between the start-up interval and the shock produces the computed 
transmission coefficient. The standard deviation of the individual responses 
serves as an error estimate. 

The dependence upon incidence angle of the acoustic transmission 

coefficient for A: * 0.001 and M = 3 waves is displayed in Fig. 3. As is 

1 s 

discussed in [19], linear theory is quite reliable at angles below, say, 
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45°. Figure 3 contains results from both spectral and finite difference 
calculations. The finite difference results were obtained with the same 
second-order MacCormack's method that was described in [19] except that 
periodic boundary conditions (rather than stretching) were employed in the y- 
direction. The finite difference grid was 64 * 16 and these calculations used 
a CFL number of 0.70. The spectral grid was 32 * 8, and the CFL number was 
0.50. Figure 3 shows that both methods produce the same results. A head-to- 
head comparison of both methods for the 0 X - 10° case is provided in Table 
III. The "exact" value is taken from linear theory [22]. Since the amplitude 
of the incident acoustic wave is so small, it should come as no surprise that 
four points in the y-direction suffice for the spectral calculation. Note 
that the standard deviations are substantially smaller for the spectral 
method. These results suggest that the spectral method requires only half as 
many grid points in each coordinate direction. 


Table III. Grid Dependence of Acoustic Transmission Coefficient 


Grid 

Finite 

Difference 

Chebyshev- 
Fourier Spectral 

16 * 4 

6.403 ± 2.652 

7.257 ± 0.587 

16 x 8 

6.427 ± 2.626 

7.257 ± 0.587 

32 x 4 

7.105 ± 0.453 

7.158 ± 0.022 

32 x 8 

7.134 ± 0.471 

7.158 ± 0.022 

32 x 16 

7.139 ± 0.497 

7.158 ± 0.022 

64 x 16 

7.163 ± 0.078 

7.157 ± 0.017 

128 x 16 

7.152 ± 0.022 


"exact" 

7.156 

7.156 
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III. PARABOLIC EQUATIONS 

The nonlinear, parabolic system formed by the incompressible, Navier- 
Stokes equations was the focus of much of the early development and 
application of spectral methods to large-scale fluid dynamical problems. 
Fourier spectral methods have been the obvious choice for the simulation of 
homogeneous, isotropic turbulence [2 A ] . For shear flows, however, non- 
periodic boundary conditions are required. So far, Chebyshev spectral methods 
have been favored for these applications [25-27]. Nevertheless, Legendre 
spectral methods are a viable alternative and of late they have been 
attracting some attention. This section will present a discussion of the 
implementation of Legendre spectral methods and will then compare them with 
Chebyshev spectral methods for the one-dimensional heat equation. This 
section will close with a description of a promising semi— implicit time- 
stepping scheme for the Navier-Stokes equations. 

Basic Legendre Spectral Concepts 

A Legendre spectral method on [—1,1] makes uses of the interpolating 
function 

N . 

u(x) = l U P (x), (47) 

n=0 

where P n (x) *-he Legendre polynomial of degree n. Closed form expres- 

sions for these polynomials are well-known, albeit clumsy. The computational- 
ly preferred way to evaluate the polynomials is through the recursion relation 

P 0 (x) = 1 
P 1 (x) = X 

and for n > 2 

n P n (x) = (2n-l)xP n _ 1 (x) - (n-1) P n _ 2 (x). (48) 

I 
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Unlike the case with Fourier and Chebyshev collocation methods, there is 
no tidy expression for the Legendre collocation points. Appeal must be made 
to the theory of numerical quadrature [28]. The presence of boundary 
conditions at both endpoints makes it desirable to include -1 and +1 in 
the set of collocation points. Subject to this constraint, the most accurate 
quadrature formula for the discrete Legendre coefficients is the Gauss-Lobatto 
rule 


u 

n 


c 

n 


N 




J-o 


P.V Uy 


n 0, 1, # **,N 


(49) 


where x Q = +1, x N = -1 and Xj for 1 < j < N-l are the roots of 
P' (x) . The weights are 


w. 


N (N+l ) P N (Xj) 


(50) 


and 



n = 0, !,*•♦, N-l 


n - N 


(51) 


The interior collocation points must be determined numerically. This 
quadrature rule yields the exact Legendre coefficients if u(x) is any 
polynomial of degree less than N. Its inverse relation is 


u, 


N ~ 

I 

n=0 


u P (*,)• 

n n j 


(52) 
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The analytic derivative of the Interpolating function in Eq. (47) is 


3u 

3x 


N 

l 


n*0 



P (x), 
n 


(53) 


where 





0 



= 0 


2n+l n 2n+5 n+2 


+ u 


n+1 


(54) 


n = N-l,N-2, ••*,0. 


Since fast transform techniques are not available for the Legendre basis 
functions, there is no particular advantage to computing 3u/3x|j by applying 
Eqs. (49), (54) and (53) rather than by following Eq. (49) with 





N 

l 


n*=0 


u 

n 


P'(X.) 
n j 


(55) 


In fact, for a collocation method it is faster still to perform this entire 
process by a single matrix-vector multiplication. For that matter the 
Chebyshev collocation differentiation operator may also be represented by a 
matrix. Timing studies [29] on the CDC Cyber 175 have indicated that even for 
N = 16, the Chebyshev matrix-multiply differentiation procedure is 

substantially faster than one based on assembly language fast transforms. 
Moreover, the matrix-multiply procedure does not suffer the sort of speed 
degradation that afflicts the transform procedure whenever N is not an 

integral power of 2. 
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The heat equation 


3u 

at 


a 2 
3 u 


ax" 


(56) 


is the natural parabolic linear model problem* The spatial domain is [-1,1] » 
the initial condition is 

u(x,0) ■» sin fix (57) 


and the boundary conditions are 


u(-l,t) - 0 

u (+1 , t ) “ 0. 


The exact solution is then 

^2 

u(x,0) = e sin ^x. 


(58) 


(59) 


The time differencing is again the classical fourth-order Runge-Kutta scheme. 

In addition to spectral collocation and series truncation solutions, we 
will also present spectral tau results. Let u n (t) ^ or n**0,l,***,N 
denote the Legendre coefficients of the tau approximation to u(x,t). The 
semi-discrete tau equations are 


du 
n 

dt 



n - 0, N-2 (60) 


with 


N 

l 

n-0 

n even 


n 


N 

* “■ 
n»l 

n odd 


0 . 


( 61 ) 


t 
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The Legendre coefficients of the approximation to the second spatial 
derivative u£ 2) (t) can be obtained from u n (t) by two applications of the 
recursion relation in Eq. (54). In this tau approximation the dynamical 
equations for the two highest-order coefficients are dropped in favor of the 
equations for the boundary conditions. Equation (61) follows from the 

property 

P (±1) = (±l) n . (62) 

n 

Since the Chebyshev polynomials also satisfy Eq. (62), the Chebyshev tau 
equations for Eq. (56) are the same as Eqs. (60) and (61). Of course, Eq. 
(29) is invoked for the derivative coefficients instead of Eq. (54). 

The results at t = 1 are given in Tables IV and V. The maximum errors 

2 

shown there have been boosted up by the factor e^ so that they represent 
relative errors. On the whole the collocation results are the best. 

Moreover, except for the truncated series results, the Legendre approximations 
are superior to the Chebyshev ones. Lanczos [30] has discussed some 
circumstances under which Legendre approximations are superior to Chebyshev 
ones. It goes almost without saying that finite difference results are far 
inferior to any of these spectral approximations. 
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Table IV. Max-timim Error for Legendre Approximations to the 


Heat Equation 


N 

Truncated Series 

Tau 

Collocation 

8 

6.65 

(-4) 

6.85 

(-4) 

2.40 

(-5) 

10 

1.72 

(-5) 

1.07 

(-5) 

1.50 

(-7) 

12 

3.06 

(-7) 

1.54 

(-7) 

1.38 

(-9) 

14 

3.50 

(-9) 

1.86 

(-9) 

4.81 

(-10) 

16 

3.88 

(-11) 

1. 15 

(-10) 

9.98 

(-11) 



Table V. Maximum Error for Chebyshev Approximations to 
Heat Equation 

the 

N 

Truncated Series 

Tau 

Collocation 

8 

2.44 (-4) 

1.61 

(-3) 

4.58 (-4) 

10 

5.76 (-6) 

2.12 

(-5) 

8.25 (-6) 

12 

9.42 (-8) 

3.19 

(-7) 

1.01 (-7) 

14 

1.14 (-9) 

3.35 

(-9) 

1.10 (-9) 

16 

1.05 (-11) 

8.39 

(-11) 

2.09 (-11) 


The time-step restriction for explicit Legendre or Chebyshev methods for 
the heat equation is very severe, scaling as 1/N^. This can pose quite a 
barrier to large-scale calculations for which a relative accuracy of 0.1% or 
so will suffice. Fortunately, many large-scale calculations can be split into 
one-dimensional, inhomogeneous counterparts of Eq. (56) and efficient implicit 
schemes are available for this linear, constant coefficient equation. They 
rely on reducing the Legendre (or Chebyshev) tau equations to a system which 
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is nearly tridiagonal* The Legendre tau equations for a Crank - Nicolson 

temporal discretization of Eq. (56) are 

. 2 A e , « ^6 i/ 

A — . r, n+2__i — . n+ 4 — 

(2n-l)(2n-'3T U n-2 1 (2n-l) (2n+3) J n (2n+3)(2n+5) n+2 

1 j _ 2e n+2 T + e n+4 T 

(2n-l) (2n-3) n-2 (2n-l)(2n+3) n (2n+3)(2n+5) n+2 

n * 2,3, •• • ,N, (63) 


where A = -At/2 with At the time-step, the coefficients on the left- 

hand side are at t + At, 


f 

n 


u R (t) + j At uf>(t), 


(64) 


and 


! 1 0 < n < N 

0 n > N 


(65) 


Equation (63) for even n plus the first of Eqs. (61) form a linear system 
which is tridiagonal except for the boundary condition equation. This is 
cheap to invert. The odd coefficients display a similar structure. The 
Chebyshev tau version of Eq. (63) is available in [15] and [31]. 

Application to Channel Flow 

Several three-dimensional Navier-Stokes algorithms have been developed 
which incorporate the quasi-tridiagonal structure of the Chebyshev tau 
equations for the second derivative in semi-implicit schemes which treat the 
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constant coefficient diffusion term implicitly [25-27]. In practice this 
device has permitted time-steps several orders of magnitude larger than the 
explicit diffusion limit. Unfortunately, the quasi-tridiagonal structure is 
lost even for a linear, variable viscosity coefficient. An effective 
iterative scheme for this more general case has recently been proposed^ - . This 
approach will be described here in its two-dimensional setting. 

The rotation form equations for two-dimensional channel flow are 

U t ~ V(V X - Uy) + P X - (PU X ) X + (UUy) y 

V t + U(V X - Uy) + Py = (WV X ) X + (PVy)y (66) 

U x + Vy = 0, 

with periodic boundary conditions in x and no-slip boundary conditions at 
y = ±1. The variable P denotes the total pressure. The viscosity y is 
presumed to depend upon y. 

A useful discretization employs Fourier series in x and Chebyshev 
series in y. The pressure gradient term and the incompressibility constraint 
are best handled implicitly. So, too, are the vertical diffusion terms 
because of the fine mesh— spacing near the wall. The variable viscosity 
prevents the standard Poisson equation for the pressure from decoupling from 
the velocities in the diffusion term. The algorithm described in [26] appears 


Hlalik, M. R., High Technology Corp., Hampton, VA; Zang, T. A., College of 
William and Mary, Williamsburg, VA and Hussaini, M. Y., ICASE, NASA Langley 
Research Center, Hampton, VA., 1982, in preparation. 
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to be a good starting point. A Crank-Nicolson approach is used for the 
implicit terms and Adams-Bashforth for the remainder. After a Fourier 
transform in x, the equations for each wavenumber k have the following 
implicit structure 

u - V 2 At(pUy)y + V 2 AtikP = ••• 

V - V 2 At(jlV y )y + V 2 ^ tPy " ••• ' (67) 

iku + v y = 0. 

Fourier transformed variables are denoted by hats, the subscript y denotes a 
Chebyshev spectral derivative, and At is the time increment. 

The algorithm in [26] was devised for constant viscosity, in which case 
the Eqs. (67) can be reduced to essentially a block-tridiagonal form. This 
cannot be done in the present, more general situation. We advocate solving 
these equations iteratively after applying a finite difference pre- 

conditioning. 

The interesting physical problems have high Reynolds number, i.e., low 

viscosity. Thus the first derivative terms in Eqs. (67) predominate. The 

effective pre-conditioning of them is crucial. Four possibilities have been 
considered. The eigenvalues of pre-conditioned iterations for the model 
scalar problem u x = f with periodic boundary conditions on [0 , 2tt ] are 

given for each possibility in Table VI. The term aAx is the product of a 
wavenumber a and the grid spacing Ax. It falls in the range 
0 < | aAx | < tt . For the staggered grid case the discrete Eqs. (67) are 
modified so that the velocities and the momentum equations are defined at the 
cell faces yj = cos(iTj/N), j=0,l,***,N, whereas the pressure and the 
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continuity equation are defined at the cell centers = cos(ir ( j- V 2 )/n) , 
j=l,»««,N. Fast cosine transforms enable interpolation between cell faces 
and cell centers to be implemented efficiently. The staggered grid for the 
Navier-Stokes equations has the advantage that no artificial boundary 
condition is required for the pressure at the walls. 

The actual eigenvalues for pre-conditioned iterations of Eqs. (67) are 
displayed in Fig. 4. The model problem estimates the eigenvalue trends 
surprisingly well considering that it is just a scalar equation, has only 
first derivative terms and uses Fourier series rather than Chebyshev 
polynomials. 

The preceding results indicate that the staggered grid leads to the most 
effective treatment of the first derivative terms. The condition number of 
the pre-conditioned system is reasonably small and no resolution is lost by a 
high mode cut-off. (Although it is possible to devise a high-mode cut-off 
which avoids the small eigenvalues shown in the figures, some of the spectral 
resolution is thereby lost.) A simple and effective iterative scheme for this 
system with its complex eigenvalues is a minimum residual method. At a 
Reynolds number of 7500 each iteration reduces the residual by almost an order 
of magnitude. 

Table VII presents a comparison of the accuracy of the Chebyshev 
discretization in y. The two codes are otherwise identical. The initial 
condition consisted of Poiseuille flow plus a small amount of a linearly 
unstable eigenmode. The table compares the computed growth rate of this 
perturbation with the theoretical, linear result after 100 time-steps. 
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Table VI. Pre-conditioned Eigenvalues for a One— dimensional 
First Derivative Model Problem 


PRE-CONDITIONING 

EIGENVALUES 

Central Differences 

aAx 

sin (aAx) 


One-sided Differences 

-i(aAx/2) 

aAx/2 

e — 

s 

in( (aAx)/2)J 

High Mode Cut-off 

( aAx 
) sin(aAx) 

0 < |aAx| < ( 2 tt /3 ) 

( 0 

(2tt/ 3) < laAx! < ir 

Staggered Grid 

(aAx)/2 


sinl (aAx)/2l 



Table VII. 

Percent Error in Growth 

Rate 

N 

Finite Difference 

Spectral 

8 

4470 

3210 

16 

337 

74.5 

32 

147 

0.097 

64 

39.5 

0.071 

128 

10.0 


256 

2.4 
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IV. ELLIPTIC EQUATIONS 

Fruitful nonlinear applications of spectral methods developed the latest 
for equations of elliptic type. Unlike hyperbolic or parabolic equations, for 
which explicit schemes can often be tolerated, ellipic equations virtually 
require implicit iterative schemes in practical situations. It was only a few 
years ago that Morchoisne [32] and Orszag [33] proposed preconditioning the 
spectral collocation equations by finite difference operators. More recently 
still, effective spectral multigrid iterative methods have been developed 
[34,35] and applied to the nonlinear potential flow problem of fluid dynamics 
[29] . These developments will be described in this section. 


Poisson's Equation 

As usual the discussion will begin with a linear model problem, but this 
time in two spatial dimensions. That problem is the Poisson's equation 




f 


( 68 ) 


on the square [-1,1] x [-1,1] with homogeneous Dirichlet 
conditions. The choice 

2 

f(x,y) = -2rr sin irx sin iry 


boundary 

(69) 


corresponds to the analytical solution 

u(x,y) = sin irx sin Try. (70) 

Both Chebyshev and Legendre spectral methods are appropriate for this 
problem. Direct solution schemes for the Chebyshev tau method have been 
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described in [31]. The same schemes also work for the Legendre tau method 
with straightforward modifications. They are basically of an alternating 
direction implicit (ADI) nature and rely on the quasi-tridiagonal form of the 
constant coefficient, one-dimensional problem. Haidvogel and Zang [31] report 
comparisons of the Chebyshev tau method with finite difference methods on 
numerous problems. They discuss both computational efficiency and accuracy. 

These direct solution schemes cannot feasibly be extended to spectral 
collocation methods because the collocation equations for the one-dimensional 
components cannot be represented by sparse matrices. However, an ADI 
iterative scheme based on finite difference preconditioning is an efficient 
method for obtaining an approximate solution. The description of this scheme 
in its general nonlinear setting begins by writing the spectral collocation 
equations as 

M(U) - 0. (71) 

Define the Jacobian 

J (U) = w < u >* (72) 

In many cases the Jacobian can be split into the sum of two operators 

J X (U) and J (U), each involving derivatives in only the one coordinate 

direction indicated by the subscript. The most straightforward ADI method is 

[al - J (V) ] [al - J (V) ]AV = aM(V), (73) 

x y 

with the approximate solution V updated by 


V + V + wAV. 


(74) 
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This is just the Douglas-Gunn version of ADI [36] . The term approximate 

factorization is commonly used for this type of scheme for the nonlinear 

potential flow problem [37]. This particular scheme is referred to as AF1 . 

For second-order spatial discretizations the term [al - J X (V)] leads to a 

set of tridiagonal systems, one for each value of y. The second left-hand 

side factor produces another set of tridiagonal systems. For spectral 

discretizations, however, these systems are full; hence, Eq. (73) is still 

relatively expensive to invert. A compromise is to replace J x and Jy with 

their second-order finite difference analogs, denoted by H and H , 

x y 

respectively: 

[al - H (V) ] [al - H (V)]AV = aM(V). (75) 

x y 

The spectral approximate factorization scheme consists of Eqs. (74) and 
(75). The choice of the iteration parameters is discussed in [29]. 


Table VIII. Maximum Error for Chebyshev Approximations to 
Poisson's Equation 


N Truncated Series Tau Collocation 


8 

2.88 

(-4) 

2.79 

(-3) 

1.17 

(-4) 

10 

6.79 

(-6) 

5.26 

(-5) 

2.33 

(-6) 

12 

1.09 

(-7) 

8.86 

(-7) 

3.12 

(-8) 

14 

1.34 

(-9) 

1.09 

(-8) 

3.27 

(-10) 

16 

1.19 

(-11) 

9.15 

(-11) 

2.73 

(-12) 
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The results for the simple model problem are presented in Tables VIII and 
IX. The trends are the same as they were for the heat equation: the 
collocation method is more accurate than tau and Legendre polynomials are more 
accurate than Chebyshev. (Since it is not practical to design a spectral 
method for PDE's using truncated series, those results have been ignored in 
this comparison.) 


Table IX. Maximum Error for Legendre Approximations to 
Poisson's Equation 


N 

Truncated Series 


Tau 

Collocation 

8 

6.04 

(-4) 

1.55 

(-3) 

1.77 

(-5) 

10 

1.69 

(-5) 

3.40 

(-5) 

2.48 

(-7) 

12 

3.05 

(-7) 

6.05 

(-7) 

2.27 

(-9) 

14 

3.82 

(-9) 

6.98 

(-9) 

1.99 

(-11) 

16 

3.85 

(-11) 

6.37 

(-11) 

3.06 

(-10) 


Spectral Multigrid Methods 

Iterative schemes for spectral collocation equations, such as AF1 , can be 
accelerated dramatically by applying multigrid concepts. This technique has 
been extensively developed for finite difference and finite element 
discretizations [39] and has recently been applied to spectral discretizations 
[34,35,29]. Briefly put, multigrid methods take advantage of a property 
shared by a wide variety of relaxation schemes - potential efficient reduction 
of the high-frequency error components but unavoidable slow reduction of the 
low-frequency components. 
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The fundamentals of spectral multigrld are perhaps easiest to grasp for 


the simple model problem 


d 2 u 


dx 


(76) 


on [0,2ir] with periodic boundary conditions. The Fourier approximation to 
the left-hand side of Eq. (76) at the collocation points is 


N/2-1 

l 


p= -N/2+1 


2 - 

P Up e J 


(77) 


The spectral approximation to Eq. (76) may be expressed as 


LU = F, 


(78) 


where 

U = ( u 0’ u l’*’ # ,U N-1^’ 


(79) 


= ( f Q’ f 1 ’ * * * » f N-l^ ' 


(80) 


and 


2 2 

L represents the Fourier spectral approximation to - d /dx . 
A Richardson's iterative scheme for solving Eq. (78) is 


V <- V + ui(F - LV), 


(81) 


where w is a relaxation parameter. On the right side of the replacement 
symbol («•) V represents the current approximation to U, and on the left it 
represents the updated approximation, The eigenfunctions of L are 
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^(p) = e 2lTiJp/N , (82) 

with the corresponding eigenvalues 

A(p) = p 2 , (83) 


where j = and p = - N/2+1, • •• ,N/2-l. The index p has a 

natural interpretation as the frequency of the eigenfunction* 

The error at any stage of the iterative process is V - U; it can be 
resolved into an expansion in the eigenvectors of L* Each iteration reduces 
the p' th error component to v(X ) times its previous value, where 

v(X) = 1 - u)X. (84) 


The optimal choice of to results from minimizing |v(X)| for 

X e [ X min> X max] » where X min 1 and X max = n2 / 4 * (° ne need not worr y about 
the p = 0 eigenfunction since it corresponds to the mean level of the 

solution, which is at one's disposal for this problem*) The optimal 
relaxation parameter for this single-grid procedure is 


SG X + X 

max min 

It produces the spectral radius 


(85) 


P 


SG 


max min 

X + X * 
max min 


( 86 ) 



36 


Unfortunately, P gG = 1 - 8/N 2 , which Implies that 0(N 2 ) iterations are 
required to achieve convergence. 

This slow convergence is the outcome of balancing the damping of the 
lowest-f requency eigenfunction with that of the highest-f requency one in the 
minimax problem described after Eq. (84). The multigrid approach takes 
advantage of the fact that the low-frequency modes (|p| < N/4) can be 
represented just as well on coarser grids. It settles for balancing the 
middle-frequency eigenfunction (|p| “ N/4) with the highest-f requency one 

(|p| =* N/2), and hence damps effectively only those modes which cannot be 
resolved on coarser grids. In Eqs. (85) and (86), ^ m ± n is replaced with 
X mid = X ( N /4)* The optimal relaxation parameter in this context is 


0) 


MG 


X + 
max 


mid 


The multigrid smoothing factor 


(87) 


_ _max mid (88) 

MG X + X 

max mid 

measures the damping rate of the high-frequency modes. In this example 
P MG = 0.60, independent of N. The price of this effective damping of the 
high-frequency errors is that the low-frequency errors are hardly damped at 
all. Table X compares the single-grid and multigrid damping factors for N = 
64. However, on a grid with N/2 collocation points, the modes for 
Ip I e [n/ 8, N/4] are now the high-frequency ones. They get damped on this 
grid. Still coarser grids can be used until relaxations are so cheap that one 
can afford to damp all the remaining modes, or even to solve the discrete 
equations exactly. For the case illustrated in Table X the high-frequency 
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error reduction in the multigrid context is roughly 250 times as fast as the 
single-grid reduction for N = 64. 

Let us consider just the interplay between two grids. A general, 
nonlinear fine-grid problem can be written 

L f (U f ) = F f . ( 89 ) 

The shift to the coarse grid occurs after the fine-grid approximation V f has 
been sufficiently smoothed by the relaxation process, i.e., after the high 
frequency content of the error V f - U f has been sufficiently reduced. The 
related coarse-grid problem is 

L C (U C ) = F c , (90) 

where 

F C = R[F f - L f (V f )] + L C (RV f ) . (91) 

The restriction operator R interpolates a function from the fine grid to the 
coarse grid. The coarse-grid operator and solution are denoted by L c and 
U c , respectively. After an adequate approximation V c to the coarse-grid 
problem has been obtained, the fine— grid approximation is corrected via 

V f + V f + P(V C - RV f ). ( 92 ) 

The prolongation operator P interpolates a function from the coarse grid to 
the fine grid. 



38 


A complete multigrid algorithm requires specific choices of the 
interpolation operators, the coarse-grid operators, and the relaxation 
schemes. These issues are discussed at length in [34,35,29] for both Fourier 
and Chebyshev multigrid methods. Numerous linear, variable coefficient 
examples are also provided there. The more interesting nonlinear examples 
from [29] are the subject of the remainder of this paper. 


Table X. Damping Factors for N ■ 64 


p 

Single-Grid 

Multigrid 

1 

.9980 

.9984 

2 

.9922 

.9938 

4 

.9688 

.9750 

8 

.8751 

.9000 

12 

.7190 

.7750 

16 

.5005 

.6000 

20 

.2195 

.3750 

24 

.1239 

.1000 

28 

.5298 

.2250 

32 

.9980 

.6000 

Application to Two-Dimensional 

Potential Flow 



Until the recent work of Streett [38], the discretization procedures for 
the potential equation were invariably based on low-order finite difference or 
finite element methods. Streett used a spectral discretization of the full 
potential equation and obtained its solution by a single-grid iterative tech- 
nique. The application of spectral multigrid techniques by Streett, et al. 
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[29] produced a dramatic acceleration of the iterative scheme. Even in its 
relatively primitive state the spectral multigrid scheme is competitive, and 
in some cases unequivocally more efficient, than standard finite difference 

schemes • 

After a conformal mapping from the surface of an airfoil to a circle the 
potential equation becomes 


where G is the reduced potential, R and 0 are the computational polar 
coordinates, and P is the fluid density. The reduced potential is periodic 
in 0 and it satisfies 



at R = 1, 


(94) 


G -*• 0 


as R 00 , 


(95) 


and the Kiitta condition. The density is given by the isentropic relation 

1 

p = [l -I=i-H2(q2 + q 2 - l)]*" 1 ; (96) 

the ratio of specific heats Is denoted by Y» and is the Mach number at 

infinity. The velocity components in the physical (r,0) plane are 

1 3 $ 

q r “ ~7T 3R 

1 3 £ 

q e “ rh 3© * 


( 97 ) 
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i0 

and the Jacobian between the complex physical plane (z = re ) and the 
complex computational plane (o = Re*-®) is 



Further details are provided in [38]. 

The spectral method employs a Fourier series representation in 0. 
Constant grid spacing in 0 corresponds to a convenient dense spacing in the 
physical plane at the leading and trailing edges. The domain in R (with a 
large, but finite outer cutoff) is mapped onto the standard Chebyshev domain 
[-1,1] by an analytical stretching transformation that clusters the 
collocation points near the airfoil surface. The stretching is so severe that 
the ratio of the largest-to-smallest radial intervals is typically greater 
than 1000. 

The flow past an NACA 0012 airfoil at 4° angle of attack and a freestream 
Mach number of 0.5 is a challenging subsonic and thus elliptic case. 
Nevertheless, the spectral solution on a relatively coarse grid captures all 
the essential details of the flow. The surface pressure coefficient from the 
spectral code MGAFSP [29] using 16 points in the radial (R) direction, and 32 
points in the azimuthal (0) direction is displayed in Fig. 5. The symbols 
denote the solution at the collocation points. For comparison, the result 
from the finite difference, multigrid, approximate factorization code FL036 
[40] is shown as a solid line. The grid used in the benchmark finite 
difference calculation is so fine (64 x 384 points) that the truncation error 
is well below plotting accuracy. The FL036 and MGAFSP results are identical 
to plotting accuracy. The spectral computation on this mesh yields a lift 
coefficient with truncation error less than 10“*. Spectral solutions on a 
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16 x 32 grid are thus of more than adequate resolution and accuracy for 
subsonic flows. 

In Fig. 6 are shown convergence histories from FL036, MGAFSP, and the 
finite difference, approximate factorization, single-grid code TAIR [41]. 
Meshes which yield approximately equivalent accuracy were chosen. The surface 
pressure results are the same to plotting accuracy, the lift coefficient is 
converged in the third decimal place, and the predicted drag coefficient is 
less than .001. (Actually, the spectral result is an order of magnitude more 
accurate than these limits, but the TAIR result barely meets them.) Figure 7 
demonstrates the improvement produced by the spectral multigrid scheme over 
the spectral single— grid method (AFSP). There is well over an order— of— 
magnitude gain in efficiency. 


V. A MIXED EQUATION 

The potential flow problem is much more difficult whenever the flow field 
contains both supersonic (hyperbolic) and subsonic (elliptic) regions. 
Nevertheless, the spectral raultigrid algorithm that succeeded for the subsonic 
flow case requires only a minor modification in order to succeed for the 
transonic (mixed) problem as well. 

The most expedient technique for dealing with the mixed elliptic- 
hyperbolic nature of the transonic problem is to use the artificial density 
approach of Hafez, et al. [42]. The original artificial density is 


P 


P “ M<$P 


(99) 


with 
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y = max{0,l - , (100) 

M 

where M is the local Mach number and 6p is an upwind first-order 
(undivided) difference. The spectral calculations employed a higher-order 
artificial density formula. The spectral method also required a weak 
filtering technique to deal with some high-frequency oscillations generated by 
the shock. Details are available in [38]. 

Flow Past an Airfoil 

A lifting transonic case is provided by the NACA 0012 airfoil at 

M = 0.75 and 2° angle of attack. A shock appears only on the upper surface 

00 

for these conditions and is rather strong for a potential calculation; the 
normal Mach number ahead of the shock is about 1.36. Lifting transonic cases 
are especially difficult for spectral methods since the solution will always 
have significant content in the entire frequency spectrum: the shock 

populates the highest frequencies of the grid and the lift is predominantly on 
the scale of the entire domain. An iterative scheme therefore must be able to 

damp error components across the spectrum. 

Surface pressure distributions from MGAFSP, TAIR, and FL036 are shown in 
Fig. 8. The respective computational grids are 18 x 64, 30 x 149, and 

32 x 192. The latter two are the default grids for the production finite 
difference codes. Spectral results obtained by trigonometrically 

interpolating the 18 x 64 grid results onto a much finer grid are included 
alongside the results at the collocation points. This reveals the wealth of 
detail that is provided by the rather coarse spectral grid. The shock 

predicted by TAIR is far more rounded and smeared than that of FL036, 
reflecting the coarser mesh and larger artificial viscosity used in the 
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former. The TAIR result shown is also only correct to one decimal place in 
lift as compared with a finer-grid result. Convergence histories for these 
three cases are shown in Fig. 9 along with the results for MGAFSP on a coarser 
grid (16 x 48). 

Flow Past a Circular Cylinder 

The MGAFSP code has recently been used by us for an extremely accurate 
determination of the critical freestream Mach number at which the potential 
flow past a circular cylinder first develops a supersonic region. This 
spectral calculation represents an alternative to the asymptotic series method 
employed by van Dyke and Guttmann [43] to arrive at the estimate 
Merit » .39823780 ± .00000001. 

The spectral multigrid potential code was used to determine the critical 
Mach number on several grids. On each of these grids calculations were 
performed at a half-dozen or so freestream Mach numbers. For each case the 
maximum local Mach number was determined from the computed solution. Then an 
extrapolation procedure was employed to ascertain what freestream Mach number 
produced a maximum local Mach number of unity. This value was recorded as the 
critical Mach number for that particular grid. An estimate of the extra- 
polation error was made to ensure consistency. These results are given in 
Table XI. 

Finally, these grid-dependent calculations of the critical freestream 
Mach number were extrapolated to the limit of infinite numerical resolution. 
The best result was obtained by assuming sixth-order convergence. The final 
estimate of the critical freestream Mach number is 

M cr -£ t » .3982415 ± .0000002. The difference between this estimate and the one 
by van Dyke and Guttmann is more than an order-of-magnitude greater than the 
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estimated errors. Nevertheless, the agreement of the two estimates to better 
than one part in 10^ is remarkable in itself. 

Note that the convergence rate of the spectral multigrid potential result 
(at least sixth-order) pertains to a quantity (critical freestream Mach 
number) which requires the fundamental solution (the potential) to be first 
differentiated and then extrapolated. Moreover, the MGAFSP code is so 
efficient that. all of the requisite calculations consumed less than 20 minutes 
of CPU time on the CDC Cyber 175 and were performed on grids with no more than 
2000 points. 

A comparable calculation by existing finite difference codes would likely 
exhibit only first-order convergence. It would be far more expensive both in 
terms of CPU time and storage, surely exceeding the central memory of a 
machine such as the CDC Cyber 175. Here then is an example which firmly 
establishes the utility of spectral methods for nonlinear, multi-dimensional 
problems . 


Table XI. Grid-dependent Critical Freestream Mach Numbers 


Grid 

Merit 

Error Estimate 

14 x 32 

.398289 

.000048 

18 x 40 

.3982514 

.0000099 

22 x 48 

.3982450 

.0000035 

30 x 64 

.3982422 

.0000007 
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Fig. 1 
Fig. 2 

Fig. 3 
Fig. - 

Fig. f 

Fig. 6. 
Fig. 7 


Figure Captions 


Typical shock-fitted time-dependent flow model in the physical 

plane. 

Post-shock dependence of the pressure response to a pressure wave 
incident at 10° to a Mach 3 shock. The solid line is the linear 
theory prediction. The circles are the spectral solution. 

Dependence on incident angle of the pressure response to a 0.1% 
amplitude pressure wave incident on a Mach 3 shock. The solid line 
is the linear theory result. Circles are spectral solutions; 

squares are finite difference solutions. 

4. Eigenvalues of the pre-conditioned matrices for semi-implicit 

channel flow when the streamwise wave number k = 5. The grid is 
32 x 17, the Reynolds number is 7500 and the CFL number is 0.10. 
Note the different scale used for the central differences pre- 

conditioning results. 

i. Spectral (triangles) and finite difference (solid line) surface 

pressures for a subcritical flow. 

. Maximum residual versus machine time for a subsonic flow. 

. Error in lift versus machine time for a subsonic flow from single- 
grid (AFSP) and multigrid (MGAFSP) spectral schemes. 
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Fig. 8. Surface pressures for a transonic flow. 

Fig. 9. Maximum residual versus machine time for a transonic flow. 
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